
from numpy import *
from ctypes import *

from ocemopy.grid import CGrid

class Gridgen(CGrid):
    """docstring for Gridgen"""
    
    _libgridgen = pydll.LoadLibrary("libgridgen.dylib")

    _libgridgen.gridgen_generategrid2.restype = c_void_p
    _libgridgen.gridnodes_getx.restype = POINTER(POINTER(c_double))
    _libgridgen.gridnodes_gety.restype = POINTER(POINTER(c_double))
    _libgridgen.gridnodes_getnce1.restype = c_int
    _libgridgen.gridnodes_getnce2.restype = c_int
    _libgridgen.gridmap_build.restype = c_void_p
    
    def generate_grid(self):
        
        # if self._gn is not None:
        #     self._libgridgen.gridnodes_destroy(self._gn)
        
        nbry = len(self.xbry)
        
        nsigmas = c_int(0)
        sigmas = c_void_p(0)
        nrect = c_int(0)
        xrect =  c_void_p(0)
        yrect = c_void_p(0)
        
        self._gn = self._libgridgen.gridgen_generategrid2( c_int(nbry), 
             (c_double * nbry)(*self.xbry), (c_double * nbry)(*self.ybry), 
             (c_double * nbry)(*self.beta),
             c_int(self.ul_idx), c_int(self.nx), c_int(self.ny), 
             c_int(0), POINTER(c_double)(), POINTER(c_double)(),
             c_int(self.nnodes), c_int(self.newton), c_double(self.precision),
             c_int(self.checksimplepoly), c_int(self.thin), c_int(self.nppe),
             c_int(self.verbose),
             byref(nsigmas), byref(sigmas), byref(nrect),
             byref(xrect), byref(yrect) )
        
        x = self._libgridgen.gridnodes_getx(self._gn)        
        x = asarray([x[j][i] for j in range(self.ny) for i in range(self.nx)])
        x.shape = (self.ny, self.nx)
        
        y = self._libgridgen.gridnodes_gety(self._gn)        
        y = asarray([y[j][i] for j in range(self.ny) for i in range(self.nx)])
        y.shape = (self.ny, self.nx)
        
        if any(isnan(x)) or any(isnan(y)):
            x = ma.masked_where(isnan(x), x)
            y = ma.masked_where(isnan(y), y)
        
        if self.proj is not None:
            if isinstance(xp, ma.MaskedArray): xp=xp.filled(nan)
            if isinstance(yp, ma.MaskedArray): yp=yp.filled(nan)
            lon, lat = proj(x, y, inverse=True)
            lon = ma.masked_where(isnan(lon), lon)
            lat = ma.masked_where(isnan(lat), lat)
            super(Gridgen, self).__init__(lon, lat, proj=self.proj)
        else:
            super(Gridgen, self).__init__(x, y)
        
    
    def __init__(self, xbry, ybry, beta, shape, ul_idx=0, \
                 focus=None, proj=None, \
                 nnodes=14, precision=1.0e-12, nppe=3, \
                 newton=True, thin=True, checksimplepoly=True, verbose=False):
        
        self.xbry = asarray(xbry)
        self.ybry = asarray(ybry)
        self.beta = asarray(beta)
        assert self.beta.sum() == 4.0, 'sum of beta must be 4.0'
        self.shape = shape
        self.ny = shape[0]
        self.nx = shape[1]
        self.ul_idx = ul_idx
        self.focus = focus
        self.nnodes = nnodes
        self.precision = precision
        self.nppe = nppe
        self.newton = newton
        self.thin = thin
        self.checksimplepoly = checksimplepoly
        self.verbose = verbose
        
        self.proj = proj
        if self.proj is not None:
            xbry, ybry = proj(xbry, ybry)
        
        self._gn = None
        self.generate_grid()
    
    def __del__(self):
        """delete gridnode object upon deletion (still memory leaks..)"""
        self._libgridgen.gridnodes_destroy(self._gn)



if __name__ == '__main__':
    from matplotlib.pyplot import *
    
    xbry = (0, 0, 10, 11)
    ybry = (11, 0, 1, 10)
    beta = ( 1, 1, 1,  1)
    
    grd = Gridgen(xbry, ybry, beta, (50, 50), ul_idx=2)
    
    plot(grd.x, grd.y, '-k')
    plot(grd.x.T, grd.y.T, '-k')
    
    show()